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THREE-DIMENSIONAL ELASTIC STRESS AND DISPLACEMENT ANALYSIS 
OF FINITE GEOMETRY SOLIDS CONTAINING CRACKS 

Jonathan Krxtig, John Gyeuenycsi, and Alexander Mendelson 
Lewis Research Center and 
U.S. Array Air Mobility R&D Laboratory 
Cleveland, Ohio 44135 


ABSTRACT . The line method of analysis is applied to the Navier- 
Cauchy equations of elastic equilibrium to calculate the displacement fields 
in finite geometry bars containing central, surface, and double-edge cracks 
under extensionally applied uniform loading. The application of this method 
to these equations leads to coupled sets of simultaneous ordinary differen- 
tial equations whose solutions are obtained along sets of lines in a dis- 
cretized region. Normal stresses and the stress intensity factor variation 
along the crack periphery are calculated using the obtained displacement 
field. The reported results demonstrate the usefulness of this method in 
calculating stress intensity factors for commonly encountered crack geom- 
etries in finite solids. 

INTRODUCTION . The main goal of fracture mechanics is the prediction 
of the load at which a structure weakened by a crack will fail. Knowledge 
of the stress and displacement distributions near the crack tip is of funda- 
mental importance in evaluating this load at failure. During the early 
development of crack mechanics most of the effort was focused on through- 
thickness cracks which could be characterized as two-dimensional. However, 
part- through cracks are the most common type of crack defect found in actual 
service conditions (ref. 1). 

Because of the geometric singularity associated with any crack type 
problem, only limited analytical work has been done in the past on these 
problems. Early theoretical solutions for three-dimensional flaw config- 
urations usually involved the discussion of cracks in infinite or semi- 
infinite solids (refs. 2 to 8) . For this reason, results for finite geom- 
etry stress intensity factors are usually given in terms of magnification 
factors applied to some convenient reference solution. In addition, con- 
siderable scatter exists in the reported results as obtained by different 
investigators (ref. 9). In our work these difficulties are avoided by 
solving the finite dimensional problems directly. 

Recently, approximate solutions of the finite geometry surface crack 
problem were obtained by the boundary integral equation method (ref. 10) 
and the finite element method (ref. 11). An alternate semi-analytical 
method suitable for the elastic solution of crack problems is the line 
method of analysis. Successful application of this method to finite 
geometry solids containing cracks has been demonstrated by Gyekenyesi 
and Mendelson (ref. 12). Although the concept of the line method for 
solving partial differential equations is not new (ref. 13), its appli- 
cation in the past has been limited to simple examples. The basis of 
this technique is the substitution of finite differences for the deriva- 
tives with respect to all the independent variables except one for which 
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the derivatives t are retained. This approach replaces a given partial dif- 
ferential equation with a system of simultaneous ordinary differential equa- 
tions whose solutions can then be obtained in closed form. These equations 
describe the dependent variable along lines which are parallel to the co- 
ordinate in whose direction the derivatives were retained. Application of 
the line method is most useful when the resulting ordinary differential 
equations are linear and have constant coefficients. 

An inherent advantage of the line method over other numerical methods 
is that good results are obtained from the use of relatively coarse grids. 
This use of a coarse grid is permissible because .parts of the solutions 
are obtained in terms of continuous functions. Additional accuracy in 
normal stress distributions is derived from the fact that they are ex- 
pressed as first-order derivatives of the displacements and these deriv- 
atives can be analytically evaluated. Inherently inaccurate numerical 
differentiation is required only for evaluating the shear stresses, but 
this presents no important loss of accuracy since they are an order of 
magnitude smaller than the normal stresses. For problems with geometric 
singularities, additional accuracy is derived from using a displacement 
formulation since the resulting deformations are not singular. 

It is the purpose of this report to present a simple and systematic 
approach to the elastic analysis of three-dimensional, finite geometry 
solids containing traction-free cracks. The need for these specific solu- 
tions has existed for a number of years in fracture toughness testing. 

REDUCTION OF THE NAVIER-CAUCHY EQUATIONS TO SYSTEMS OR ORDINARY 
DIFFERENTIAL EQUATIONS . Within the framework of linearized elasticity 
theory, the equations of elastic equilibrium in terms of displacements 
are 
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For a finite geometry solid with rectangular boundaries, we construct 
three sets of parallel lines (fig. 1(a)). Each set of lines is parallel 
to one of the coordinate axes and thus perpendicular to the corresponding 
coordinate plane. An approximate solution of equation (1) can be ob- 
tained by developing solutions of ordinary differential equations along 
the x-dlrectional lines. As seen in the figure, t'ner are a total of 
i “ NY x NZ such lines where NY is the number of lines along the y~ 
direction and NZ is the number of lines along the z-direction in a given 
plane, respectively. We define the displacements along these lines as 
u l> u 2> * • •» u i * derivatives of the y-directional displacements 

on these lines with respect to y are defined as v '|i» v'^, • • •» v '[i» 
and the derivatives of the z-directional displacements with respect to 
z are defined as w'^, ...» w'| t . These displacements and deriva- 

tives can, then be regarded as functions of x only since they are vari- 
ables on x-directional lines. When these definitions are used, the 
ordinary differential equation along a generic line ij (a double sub- 
script is used here for simplicity of writing) in figure 1(h) may be 
written as C 
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Similar differential equations are obtained along the other x-directional 
lines. Since each equation has the terms of the displacements on the sur- 
rounding lines, these equations constitute a system of ordinary differential 
equations for the displacements u-j., U 2 » . . . , u A . 


The set of i second order differential equations represented (5) 
can be reduced to a set of 2i first order differential equations by treating 
the derivatives of the u’s as an additional set of i unknowns, ivfc. , 
defining 
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written as a single first order 


R(x) (8) 

where U and R are column matrices of 2\ elements each and A^: is a 
2.i x 2i matrix of the constant coefficients appearing in equations (5) and 
(7). 

In a similar manner, to solve equations (2) and (3) , ordinary dif- 
ferential equations are constructed along the y- and z-directional lines 
respectively. These equations are also expressed in an analogous form 
to equations (8) J they are 

^=A 2 V + S(y) (9) 


The resulting 2t equations can now be 
matrix differential equation ' 


™ = A U + 
dx 1 


~ » A 3 W + T(z) (10) 

Equations (8) to (10) are linear first-order ordinary matrix dif- 
ferential equations. They are, however, not independent, but are coupled 
through the vectors, R, S and T whose components are given by equations 
similar to (6). The elements of the coefficient matrices A^, A 2 » and 
A3 are all constants, being functions of the mesh spacing and Poisson's 
ratio only. 

Noting that a second-order differential equation can satisfy only 
a total of two boundary conditions and since three-dimensional elasticity 
problems have three boundary conditions at every point of the bounding 
surface, some of the boundary data must be incorporated into the surface 
line differential equations. Hence, conditions of normal stress and dis- 
placement are enforced through the constants of the homogeneous solutions 
while shear stress boundary data must be incorporated into the differential 
equations of the surface lines. The application of the specified shear 
conditions permits the use of central difference approximations when sur- 
face line differential equations are constructed. The details of con- 
structing these equations are found in reference 14. 

SOLUTION OF THE SYSTEMS OE ORDINARY DIFFERENTIAL EQUATIONS . The 
systems of ordinary differential equations (8) to (10) can be solved 
by any of a number of standard techniques. The method used herein was 
basically the matrizant or Peano-Baker method of integration (ref. 15). 

FqV' ^Ykfl&tpn (8) the solution can be written as 
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with similar solutions for equations (9) and (10). U(0) is the initial 

value vector, determined from the boundary conditions. The conversion of 
given boundary data into required initial values is discussed in more 
detail in reference 14. 

A v 

The matrizant el is generally evaluated by its matrix series. 

For larger values of x, when convergence becomes slow, additive formulas 
may be used. In addition, similarity transformations can be used to diag- 
onalize the matrix A^. These various techniques for improving the ac- 
curacy are discussed in detail in reference 14. 

Since equations (8) to (10) and their boundary conditions are highly ' 
coupled, it is generally impossible to directly evaluate their solutions. 
Thus, a successive approximation procedure must be employed where assumed 
values must be used initially for the required unknowns. The cyclic 
resubstitution of the obtained solutions into the coupling vectors and 
the boundary conditions will usually converge to the correct solution, 
depending mainly on the accuracy to which the required matrizant can be 
evaluated. 

Once the successive approximation procedure has converged and the 
displacement field in the body has been calculated, the normal stress 
distributions can be obtained directly by using the stress-displacement 
equations. The shear stresses, however, can be evaluated only through 
finite difference approximations for the required displacement gradients. 

STRESS INTENSITY FACTOR . The stress intensity factor Kj was at 
first obtained from the calculated stresses and displacements by extending 
the usual definition 
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to discrete data, where R Is measured from and is normal to the crack 
front and n is the singularity. It was found, however, that due to the 
coarseness of the grid used, the usual plotting and extrapolating techniques 
gave results that were erratic and of questionable accuracy. This was 
compounded by the fact that the precise crack tip location is not really 
known except that it is approximately midway between two lines, one of 
which has zero displacement specified in the crack plane and one of which 
has zero stress specified. It was found, however, that by using two terms 
in the stress and displacement series expansions around the crack tip, 
good results could be obtained even with the coarse grid used. Furthermore, 
this also permitted us to determine the actual crack tip location from the 
computed results. The method utilized is as follows. We taka 
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where ct is a function of Poisson's ratio, n was assumed to be -1/2 
and r is the crack edge position correction measured from the originally 
assumed midpoint position. Using displacement data from three adjacent 
nodes to the crack edge in equation (13), values of a K^, and r 

are calculated for each value of z, with R also measured from the 
half-way point between nodes specifying boundary stresses and displacements, 
respectively. Substituting values of Lj/Kj and r into equation (14), 
we can calculate Kj as a function of the corrected crack edge distance, 
p = R - r. A plot of In Ki versus -i/p” as "i/p" 0 can then be used to 

obtain Kj. In a similar manner, a can now be calculated from equation (13), 
where the corrected crack edge distance with the displacement data is 
p = R + r. 

APPLICATION TO TENSILE FRACTURE SPECIMENS CONTAINING CRACKS . A great 
amount of experimental work has .been done in fracture mechanics (ref. 16) 
through the use of crack-notched specimens. In the past, many different 
types of specimens have been used to determine a material's fracture toughness. 
The most common early specimens employed in these tests were the center- 
cracked and double-edge-notched bar specimens. Figures 2(a) and 3(a) show 
the finite rectangular bars with through-thickness, traction-free central 
and double-edge cracks, respectively. Because of the symmetric geometry and 
loading, only one-eighth of each bar has to be discretized as shown in 
figures 2(b) and 3(b). 

NUMERICAL RESULTS . - Center-Cracked Tensile Fracture Specimen. The 
solution of this problem was obtained by using two different sets of lines 
along the coordinate axes so that the convergence of the finite difference 
approximations could be checked. In a given direction, uniform line spacing 
was used in all computations with no other restriction being placed on the 
selection of the grid size. The crack edge location with respect to the 
imposed grid was initially assumed to be halfway between nodes specifying 
normal stress and displacement boundary conditions, respectively. Sub- 
sequently, using the obtained near crack tip stresses and displacements, 
a more accurate crack edge location was established for calculating the 
stress intensity factor. This approach was considered acceptable since 
the results from the two sets of lines at corresponding points did not 
change, although the crack edge to node distance was considerably decreased 
for the finer mesh. The successive approximation procedure required for 
decoupling the three sets of ordinary differential equations was terminated 
when the difference between successively calculated nondimens ionalized 
displacements, which are of the order of unity, at every point was less 
than a present value (10“®) . As expected, the convergence rate of this 


7 


successive approximation procedure was greatly dependent on the initial 
guess for the required unknowns in the coupling vectors and boundary condi- 
tions, For maximum computer efficiency, displacement data obtained from the 
use of coarse grids was interpolated to obtain improved starting values 
for the computations involving the final spacing of lines . The required 
initial quantities for the preliminary coarse grid calcualtions were taken to 
be zero in our work. All calculations were performed on a UNIVAC 1100/40 
computer, using double precision arithmetic. 

For the selected geometry, the crack opening displacements and normal 
stresses from our analysis and those from Raju's finite element method 
(ref. 17) are compared in Table I. Although in our solution of the gross 
displacement and stress fields, the minimum corrected crack edge distance 
is p = 0.042c, crack opening displacements and stresses can be calculated 
from equations (13) and (14) at any value of p for which these equations 
are assumed to be accurate (p ** 0.40c or less in this problem). As seen 
from Table I, there is good agreement in most displacements, with the normal 
stress at the surface showing the greatest difference. 

The dimensionless crack opening displacement is shown in figure 4. 
Agreement with the finite element results is seen to be very good. It is 
noteworthy that the results correspond to elliptical crack profiles in all 
cases . 

An indication of the accuracy of our technique for computing the stress 
intensity factor is seen in figure 5. This figure shows the stress intensity 
factor variation across the bar thickness. The results obtained in ref- 
erence 17 using a finite elements method for a geometry almost identical 
to one of the geometries in this paper is also shown. It Is seen that very 
good agreement is obtained between two completely different methods. In 
addition, Isida's plane solution (ref. 18), corrected for finite width and 
length, is also shown in this figure for comparison. Note that these results 
indicate a small increas in Kj at the surface with increasing bar thickness. 
Interestingly, for bars with t > 3c, K-j- increases gradually with z, 
reaching a maximum near z ■ 0.85t, and then decreases rapidly to its surface 
value . 

Surface Crack Tensile Fracture Specimen . Figure 6 shows a finite geom- 
etry bar containing a traction free rectangular surface crack. Because of 
the symmetric geometry and loading, only one-fourth of the bar has to be 
discretized as shown in figure 6(b). 

Selected results of the dimensionless surface crack opening displacements 
are shown in figure 7. Note that the crack opening increases rapidly with 
crack depth for 0.21 <_ a/t < 0.87, slightly exceeding even the surface 
crack displacement of a through-thickness crack at a/t = 0.87. The plane 
strain solution for a finite width center cracked bar is also shown in figure 7 
for reference. Final displacement values in this report were obtained from 
a set of 100, 140, and 140 x-, y-, and z-direetional differential equations, 
respectively. A typical computer run for this system of equations takes 
approximately 30 to 40 minutes of CPU time and 720 K bytes of storage. 
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In order to show the singularity of the stresses, the y-directional 
normal stress in the crack plane is plotted in figure £ for a/t — 0.536. 

The results clearly indicate the singular nature of o along the crack 
periphery. ^ 

Double Edge Crack Tensile Fracture Specimen . Our last example is the 
finite bar with double edge cracks . The crack opening displacements for this 
problem are presented in figure 9. The stress intensity factor variation 
as a function of bar thickness is shown in figure 10. In both, results 
from the finite element method are shown for comparison. Agreement is 
again excellent. 

CONCLUSIONS . The line method of analysis presented affords a practical 
way for analysis of three-dimensional crack problems, at least for bodies 
with reasonably regular boundaries. Because parts of the solution are ob- 
tained as continuous functions along the lines chosen, relatively good 
accuracy can be obtained with coarse grids. Results of the analysis include 
the displacements and normal stresses at every node inside the body from 
which the stress intensity factor variations were easily calculated. In 
addition it should be noted that the common semi-elliptical surface crack 
problem could also be analyzed by merely changing the boundary conditions 
at certain nodes in the crack plane. Introduction of plasticity into the 
analysis could also be accomplished by changing the coupling terms in 
equations (8) to (10). Since these have to be determined by an interative 
process in any case, it would seem possible to solve the elastoplastic prob- 
lem by a simple extension of the present method. Whether this approach is 
practical requires further investigation. 

REFERENCES 


1. Srawley, John E.J and Esgar, Jack E.: Investigation of Hydrotest 

Failure of Tniokol Chemical Corporation 260-Inch-Diameter SL-1 Motor 
Case. NASA TM X-1194, 1966. 

2. Irwin, G. R. : Crack-Extension-Force for a Part-Through Crack In a 

Plate. J. Appl. Mech. vol. 29, no. 4, Dec. 1962, pp. 651-654. 

3 t Smith, Fred. W. : Stresses Near a Semi-Circular Edge Crack. Ph.D. 

Thesis, Univ. Washington, 1966. 

4. Smith, F. W. ; and Alavi, M. J.: Stress Intensity Factors for a Part- 

Circular Surface Flaw. Proceedings of the First International Pressure 
Vessel Conference, Holland, 1969. 

5. Thresher, R. W.j and Smith, F, W,: Stress-Intensity Factors for a 

Surface Crack in a Finite Solid, J. Appl. Mech,, vol, 39, no. 1, 

Mar. 1972, pp, 195-200. 

6. Shah, R. C.; and Kobayashi, A, S.: Stress Intensity Factor for an 

Elliptical Crack Under Arbitrary Normal Loading. J. Eng. Fract. 

Mech., vol. 3, July 1971, pp. 71-96. 

7. Shah, R. C.J and Kobayashi, A. S»: On the Surface Flaw Problem, 

The Surface Crack: Physical Problem and Computational Solutions, 

Am, Soc. Mech, Eng., 1972, pp. 79-124. 


9 


8. Hartranft, R. J.; and Sih, G. C.: Alternating Method Applied to Edge 

and Surface Crack Problems. Methods of Analysis and Solutions of 
Crack Problems: Recent Developments in Fracture Mechanics; Theory 

and Methods of Solving Crack Problems. Leiden, Noordhoff International 
Publ. , 1973, pp. 179-238. 

9. Keays, R. H. : A Review of Stress Intensity Factors for Surface and 

Internal, Cracks. ARL/SM-Rept-343, Aeronautical Res. Labs., 1973. 

10. Cruse, T. A.: Numerical Evaluation of Elastic Stress Intensity Factors 

by the Boundary-Integral Equation Method. The Surface Crack: Physical 

Problems and Computational Solutions. Am. Soc. Mech. Eng., 1972, 
pp. 153-170. 

11. Marcal, P. V.: Three-Dimensional Finite Element Analysis for Fracture 

Mechanics. The Surface Crack: Physical Problems and Computational 

Solutions. Am. Soc. Mech. Eng., 1972, pp. 187-202. 

12. Gyekenyesi, J. P.; and Mendelson, A.: Three-Dimensional Elastic Stress 

and Displacement Analysis of Finite Geometry Solids Containing Cracks. 
Int. J. Fract., vol. II, no. 3, June 1975, pp. 409-429. 

13. Jones, D. J.; South, J. C.; and Klunker, E. B.: On the Numerical 

Solution of Elliptic Partial Differential Equations by the Method of 
Lines. J. Comput. Fhys., vol. 9, 1972, pp. 496-527. 

14. Gyekenyesi, J. P.: Solution of Some Mixed Boundary Value Problems of 

Three-Dimensional Elasticity by the Method of Lines. Ph.D. Thesis, 
Michigan State Univ. , 1972. 

15. Frazer, R. A.; Duncan, W. J.; and Collar, A. R. : Elementary Matrices 

and Some Applications to Dynamic and Differential Equations. Cam- 
bridge Univ. Press, 1938. 

16. Sih, G. C.: Handbook of Stress Intensity Factors. Inst. Fract. 

Solid Mech., Lehigh Univ., 1973. 

17. Raju, I. S.; and Newman, J. C. , Jr.: Three-Dimensional Finite-Element 

Analysis of Finite Thickness Fracture Specimens. NASA TN D-8414, 

1977. 

18. Isida, M. : Effect of Width and Length on Stress Intensity Factors of 

Internally Cracked Plates Under Various Boundary Conditions. Int. 

J. Fract. Mech., vol. 7, Sept. 1971, pp. 301-316. 


E-9263 


TABLE L - DIMENSIONLESS CRACK OPENING DISPLACEMENTS Ev/o Q c AND DIMEN- 
SIONLESS y-DIRECTlONAL 5TRESSES Oy/o 0 FOR A RECTANGULAR BAR UNDER 
UNIFORMTENSION CONTAINING A THROUGH-THICKNESS CENTRAL CRACK 


[ALL DATA IN y=0 PLANE WITH p - CORRECTED CRACK EDGE DISTANCE. 1 


FINITE ELEMENT METHOD - REFERENCE 10 
C*1.D; W '2.0; L = 1. 75; t “ 3. 0; v = 1/3 

LINE METHOD 
c-1.0; W-L92; L-1.69; t 

» 2. 90; V 1/3 
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_Ev 

0 0 c 
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°0 

P 

Ev 

°o c 
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°0 
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0.0132 

0.3991 

9.126 

0.0132 

0.384 

8.975 


.0264 

.5626 

6.650 

.0264 

.541 

6.552 


.040 

.6910 

5.381 

.040 

.666 

5.599 


080 

.9698 

3.852 

.080 

.937 

3.996 


120 

1.1758 

3.245 

.120 

1.142 

3.330 


.160 

1.3475 

2736 

.160 

1.314 

2912 


.200 

1.4928 

2 434 

.200 

1.463 

2602 


.400 

L9935 

1.703 

.400 

2023 

1.888 


.600 

.920 

2.3120 

1.295 

.600 

.920 

2390 

2640 

1.3a 

.532 


LOCO 

2.5347 

.531 

1.000 

2.651 


t 

0.0132 

0.4293 

7.748 

0.0132 

0.441 

8.463 

2 

.0264 

.6102 

5.259 

.0264 

.622 

6.121 

.040 

.7552 

4.046 

.040 

.766 

5.030 


.080 

1.0S01 

2821 

.080 

1.081 

3.667 


.120 

L3282 

2426 

.120 

1.321 

3.048 


.160 

1.5391 

2100 

.160 

1. 524 

2 664 


.200 

1.7197 

1.989 

.200 

1.679 

2392 


.400 

2.3433 

1.533 

.400 

2.385 

1.678 


.600 

2.7412 

L 253 

.600 

2830 

1.210 


.920 



.920 

3.140 

.656 


1,000 

3.0167 

.703 

1.000 

3.157 




(a I THREE SETS OF LINES PARALLEL TO X-, AND z- 
COORDINATES AND PERPENDICULAR TO CORRESPONDING 
COORDINATE PLANES. 



(b) SET OF INTERIOR LINES PARALLEL TO x-COORDINATE. 
Figure 1. - Sets of lines parallel to Cartesian coordinates. 




(a) RECTANGULAR BAR WITH THROUGH-THICKNESS CENTRAL 
CRACK. 



(b) DISCRETIZED REGION OF RECTANGULAR BAR WITH 
THROUGH-THICKNESS CENTRAL CRACK. 


Figure 2. - Rectangular bar with through-thickness central 
crack under uniform tension. 



(a) RECTANGULAR BAR WITH THROUGH-THICKNESS DOUBl£ 
EDGE CRACKS. 



(b) DISCRETIZED REGION OF RECTANGULAR BAR WITH DOU5H- 
EDGE CRACKS. 

Figure 3. - Rectangular bar with through-thickness double-edge 
cracks under uniform tension. 






Figure 4. - Crack opening displacement for center-cracked bar under 
uniform tension. 



Figure 5. Stress-intensity (actor variation as a (unction of bar thick- 
ness tor a ce ( -T-cracked rectangular bar under uniform tension. 






(b) DISCRETIZED REGION OF BAR WITH RECTANGULAR 
SURFACE CRACK. 



Figure 7. - Surface crack opening displacement variation as 
a function of crack depth for a rectangular bar under uni- 
form tension. 


Figure 6. - Bar with rectangular surface crack under uniform tension. 
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Figure 9. - Crack opening displacement for rectangular bar under 
uniform tension containing through-thickness double-edge 
cracks. 



z/t 


Figure 10. - Stress intensity factor variation as function of bar thick- 
ness for rectangular bar under uniform tension containing double- 
edge cracks. 
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